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> • ABSTRACT 

o: 

^ . Fragmentation process in a cylindrical magnetized cloud is studied with the nested 

grid method. The nested grid scheme use 15 levels of grids with different spatial reso- 
lution overlaid subsequently, which enables us to trace the evolution from the molecular 
cloud density ~ lOOcm"'^ to that of the protostellar disk ~ 10^°cm~^ or more. Fluc- 
^ . tuation with small amplitude grows by the gravitational instability. It forms a disk 

O ! perpendicular to the magnetic fields which runs in the direction parallel to the major 

axis of the cloud. Matter accrets on to the disk mainly flowing along the magnetic fields 
and this makes the column density increase. The radial inflow, whose velocity is slower 
. than that perpendicular to the disk, is driven by the increase of the gravity. While the 

Q^ \ equation of state is isothermal and magnetic fields are perfectly coupled with the matter, 

"^ I which is realized in the density range of p ^ lO^'^cm"^, never stops the contraction. The 

Q^' structure of the contracting disk reaches that of a singular solution as the density and 

Q . the column density obey p{r) oc r~^ and a{r) oc r~^, respectively. The magnetic field 

^ \ strength on the mid-plane is proportional to p{rY^'^ and further that of the center (Be) 

52 I evolves as proportional to the square root of the gas density (oc p^^). It is shown that 

isothermal clouds experience "run-away" collapses. The evolution after the equation of 
state becomes hard is also discussed. 
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1. Introduction 

Molecular clouds are often observed as filaments. 
One example is p Ophiuchi cloud, in which two fil- 
aments extends from p Oph main cloud [l ~ 353°, 
h ~ 17°) and p Oph East cloud {I ~ 354°, b ~ 16°), 
respectively. The filaments have full lengths of 10- 
17.5 pc, although their widths are as narrow as 0.24 
pc (Loren 1989). In the Taurus molecular cloud (this 
cloud looks like a couple of filaments as a whole) , be- 
tween Hcilcs Cloud 2 and L1495, a filamentary dark 
cloud (B213) is seen. Its size is approximately 20' x 
3°~ (1.2 pc X 11 pc) (Fukui & Mizuno 1991). The 
axial ratio reaches ^ 10 in this case, too. Filaments 
are found even in giant molecular clouds (such as the 
Orion L1641 cloud: Bally 1989) which are known as 
sites of massive star formation. In Figure 10 of Bally 
(1989), the topology of the magnetic field is shown 
for L1641, which indicates the magnetic fields in the 
cloud are running in the direction parallel to the ma- 
jor axis of the cloud, while just outside of L1641 
toroidal magnetic fields are found lapping the fila- 
ment. 

Collapse and fragmentation in the cylindrical cloud 
are studied by Bastien and collaborators for non- 
magnetic clouds (Bastien 1983); Tomisaka (1995) and 
Nakamura, Hanawa, & Nakano (1994) for magnetized 
clouds. Summarizing result by Tomisaka (1995; here- 
after referred to as Paper I) and Nakamura et al 
(1994), the process of fragmentation and collapse in 
the cylindrical magnetized isothermal cloud is as fol- 
lows: The eigen-mode which has the most unstable 
growth rate predominates even from random fluctua- 
tions with small amplitudes and the cylindrical cloud 
fragments into prolate spheroidal condensations sep- 
arated by every A ~ 20cs/(47rG'/9c)^^^ where Cs and 
Pc represent, respectively, the isothermal sound speed 
and the gas density on the symmetric axis of the cylin- 
der. This condensation contracts mainly in the direc- 
tion parallel to the magnetic field (axis of the cylin- 
der). Finally forms a disk which runs perpendicu- 
larly to the magnetic field. The disk, which continues 
to contract, has the density distribution similar to 
that of a singular solution for isothermal hydrostatic 
sphere (Chandrasekhar 1939): psing = c^/27rGr^ in 
the r-direction. In z-direction, the density distribu- 
tion is much different from the above but the scale 
height in z-dircction is much smaller than that of r- 
direction. Accretion occurs mainly in the direction 
parallel to the magnetic field [v^ ^ Vr). 



Since as the contraction proceeds the density in 
the center of the disk increases, spatial resolution is 
needed to see the evolution completely. First method 
to resolve it is to increase the mesh number. However, 
time necessary for the calculation increases at least 
proportional to ex (mesh number in one dimension)^ 
(see Appendix of Paper I). Thus, it is difficult to use 
a grid with meshes in one dimension TV ^ 1000 even 
using massive supercomputers. Here, we adopt the 
"nested grid method" (Ruffert 1992; Yorke, Bohden- 
heimer, & Laughlin 1993). In this scheme, several 
levels to grids are prepared: a fine grid which covers 
only the disk center, and a coarse grid which covers 
the whole cloud. The grids system is overlaid with 
each other and connected as the boundary conditions 
for a fine grid is determined by a coarse grid and con- 
trarily physical quantities of a coarse grid are calcu- 
lated by those of an overlaid fine grid. This idea was 
originally developed by Berger & Oliger (1984) and 
Berger & ColeUa (1989). 

2. Models and Numerical Method 

Since for the dense gas found in interstellar clouds 
with ~ lOK, the equation of state is well approxi- 
mated with the isothermal one (Larson 1985), we as- 
sume here that the gas obeys the isothermal equa- 
tion of state. While the gas density does not ex- 
ceed n ^ 10^°cm"'^, the decouple process of the mag- 
netic fields due to Joule loss is not effective (Nakano 
1988). Further the time scale of the ambipolar dif- 
fusion/plasma drift (the time scale that the neutral 
molecules fiow toward the cloud center across the 
magnetic fields and thus charged ions) is estimated as 
^ 10 times the free-fall time (Nakano 1988). Since the 
time scale in which the gravitational instability grows 
is shorter than that of the ambipolar diffusion, it is 
reasonable if we assume the magnetic fields are per- 
fectly frozen-in the matter even for the neutral mat- 
ter. Thus, we study the evolution of the cloud under 
the assumptions of isothermal equation of state and 
ideal magnetohydrodynamics (MHD). 

As in Paper I, we assume the length of the cylindri- 
cal cloud is much longer than the width of it and we 
apply a periodic boundary condition in the direction 
of the axis of the cloud. The length of the numerical 
box is a priori arbitrary. Since it is shown, in Paper 
I, that only the wave with the maximum growth rate 
appears even from a white noise, we restrict ourselves 
to a numerical box whose size is equal to the most un- 



stable wavelength known by the linear stability theory 
(e.g., Elmegreen & Elmegreen 1978; Nagasawa 1987) 

as 

2O7C, 



A. 
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where pc means the density on the symmetric axis of 
the cylinder and 7 ~ 1 represents a numerical fac- 
tor obtained from numerical calculation (Nakamura, 
Hanawa, & Nakano 1993). 



2.2. Basic Equations 

Basic equations are the unsteady MHD equation 
and the Poisson equation for the gravitational poten- 
tial. We assume no toroidal fields (neither magnetic 
fields nor velocity). In the cylindrical coordinate (2;, 
r, 0)q| with 9/90 — 0, they are expressed as follows: 
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2.1. Initial Condition 

The initial condition is the same as Paper I. That 
is, we assume that the isothermal cylindrical cloud 
which is in a magnetohydrostatic equilibrium and im- 
mersed in an external pressure, Pext = c^sPs, where Cg 
and ps represent the isothermal sound speed in the 
cloud and the gas density on the surface of the cloud. 
To mimic the situation that fluctuations with small 
amplitudes grow with the gravitational instability, we 
add density inhomogeneity with a small amplitude. 
When a ratio of the magnetic pressure B1/^tt to the 
thermal one is assumed constant, a/2 = Bl/Sir/p = 
constantR, the density distribution in hydrostatic bal- 
ance is given by 
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■nGpcT 



2cf(l + a/2) 



(2-2) 



where no toroidal component of the magnetic fields 
exists, Brf, = 0. Scaling in a unit such as AttG = Cs = 
Pcxt = 1, this is rewritten as 



Po{r) = Pel 



8 H-a/2 



(2-3) 



where the distance is normalized with H = Cs/(47rGps)^^^, 
and the density is normalized with that on the surface 
of the cloud ps — Pcxt/c^. The radial density distri- 
bution is seen in Figure 1 of Paper I. Parameters used 
here are summarized in TablqlJ. 

As a fluctuation added to pq, we assume a sinu- 
soidal function with a wavelength l^ , 

p{z, r) = pQ{r){l - 6) cos{2ttz/1,), (2-4) 

where the relative amplitude S is taken small enough 
(5=10-2). 



^a^'^ is the parameter which represents the ratio of the Alfven 
speed and the isothermal sound speed. 
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where the variables have their ordinary meanings. 
Equation (2-5) is the continuity equation; equations 
(2-6) & (2-7) arc the equations of motion. The in- 
duction equations for the poloidal magnetic fields are 
equations (2-8) & ( ^-91 ). The last equation (2-lC) is 
the Poisson equation. 



'^(z,r,<p) is the right-handed coordinate system. Thus, the 2- 
axis is illustrated as a horizontal axis in figures contrary to a 



common usage. 
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a 


F 


5 


N 


Iz 


Eq. of State 


A 


1 


100 


10-2 


14 


1.935 


Isothermal 


B 


0.1 


100 


10-2 


14 


2.21 


Isothermal 


C 


1 


100 


10-2 


14 


1.935 


Composite'^ 


D 


1 


100 


10-2 


14 


1.935 


Composite*^ 



Table 1: Model Parameters. 



"Isothermal for p < 10* ps and p oc p*' ^ for higher density. 
''The same as Model C but the critical density Pcrit is taken as 
106p,. 



2.3. Numerical Scheme 

As is seen in Paper I, the density distribution which 
is finally achieved in the gravitationally contracting 
disk is rather singular as p{r) ~ Ar^^. Therefore, 
near the center of the disk, the density increases 
monotonically and thus the density scale-height de- 
creases. Any numerical schemes do not guarantee 
the stabilities after the quantities differ much in adja- 
cent cell-to-cell, although a specific value of the ratio 
over which the calculation becomes unstable depends 
upon a scheme. In Paper I, a nonuniformly spaced 
grid was used in which the grid spacing in the z- 
direction increases as the distance from the disk mid- 
plane increases was used. However, this leads to a 
mesh whose shape is far from the square (Az = Ar), 
i.e., Az ^ Ar (near the pole r — but far from the 
disk center) or Az ^ Ar (near the disk midplane but 
far from the z-axis). Unnecessarily close spacing in 
the above has disadvantages since the numerical time 
step is severely limited there. Further, it is thought 
that schemes using the square grids are more robust 
than those using rectangle grids (Az j^ Ar). 

In the nested grid scheme, all cells are square. In 
this scheme, several levels of grid systems are pre- 
pared; a fine grid covers the disk center and a coarse 
grid is for a whole cloud. Respective grids are over- 
laid as a co-centered fashion with each other as Fig- 
ure 1. Fifteen levels of grids are used in the present 
paper, which are named as LO (the coarsest), LI, 
L2, . . ., L14 (the finest). Since the mesh spacing 
Az = Ar = A/i of the n-th level is just a half of 
that of the (n — l)-th level, the spacing of the L14 
grid is equal to 1/21"* = 1/16364 of that of LO, that 
is, Ahi4 = Aft,o/16364. Thus the scheme has a wide 
dynamic range in the space dimension. Since we take 
64 meshes in one-dimension for each levels (Ln for 
n = 0, 1, . . . , 14), the grid spacing in L14 corresponds 
to 1/1048576 ~ 1/10^ of the size of the coarse grid 
LO ih). 

Boundary values for a fine grid (Ln) are determined 
by interpolation of the values of a coarser grid (Ln — 1 ) 
and the physical quantities in a coarse cell (Ln— 1) are 
calculated by averaging those in four fine cells (Ln) 
overlaid onto the coarse mesh. Detailed procedure 
is written in Appendix. The simulation simulation 
begins with five levels (LO, LI, L2, L3 & L4). Depth 
of the grid levels is increased when spatial resolution 
is necessary, upto L14. 

Numerical methods applied to each levels of grids 



are the same as those of Paper I. (In other words, 
the scheme adopted in Paper I is identical with the 
nested grid scheme consisting of LO only.) Unsteady 
MHD equations are solved using van Leer's mono- 
tonic interpolation (1977) and the constrained trans- 
port method by Evans & Hawley (1988). We adopt 
MILUCGS (modified incomplete LU decomposition 
preconditioned conjugate gradient squared method: 
Meijerink & van der Vorst 1977; Gustafsson 1978). 
We compared the numerical result by the nested grid 
scheme with that of Paper I. While the maximum den- 
sity is less than ~ W^-^ps, the density distributions 
are identical except for regions near the boundaries. 
After the maximum density exceeds ~ 10^ ps, numer- 
ical oscillation appeared near the center of the disk in 
Paper I. However, the nested grid can resolve smaller 
structures than Paper I seen in the nest section. 

3. Results 

Model A in the present paper is identical to Model 
B of Paper I. Figure 2 shows the structure of the cloud 
at i = 1.367 when pc — 10^ ps- This corresponds to 
Figure 3 of Paper I, which is the final state that can 
be reached by a scheme using 400x400 meshes with- 
out the nested grid technique. Gas falls down along 
the magnetic field and forms a disk running perpen- 
dicularly to the magnetic fields (Fig. 2 a). Magnetic 
field lines are squeezed by a dense gas disk and form 
a valley in the center. Figure 2b shows the structure 
seen using the L2 grid, whose spatial resolution is 4 
times finer than LO. The disk whose shape is con- 
caved is captured well in L2. Although the center of 
the disk looks like a sphere in LO, this is due to a 
low resolution of LO. The maximum density attained 
is also restricted in LO as ^ lO^Ps, while that in L2 
reaches ^ lO-^-^p^ (pc ~ W^Ps for Ln n > 3). The 
scheme with evenly spaced meshes adopted in Paper 
I could not proceeds further. 

In Figure 3, we compare the structure at i = 1.372 
when Pc = W^ps (a) and that at i = 1.373 when 
Pc = lO^ps (6). Since the density increases monoton- 
ically in a short time scale, it seems better to choose 
the central density pc to indicate the time evolution, 
instead of the time passing after the simulation be- 
gins. Comparing (a) and (6), it is clearly seen that the 
density distribution only near the disk center evolves; 
The gas accrets from the region near the axis fiowing 
perpendicularly to the disk; The evolution time scale 
for the disk far from the z-axis is much longer. 



Figure 4 a and b shows the disk covered in the L9 
grid. It is seen that the disk continues to contract in 
the very center. Comparing the directions of veloci- 
ties and magnetic fields, we can see the gas flows in 
the direction parallel to the magnetic field outside the 
disk (a region with p <^ 10^). This means that the gas 
flow is controlled by the magnetic fleld and the mag- 
netic field lines are not squeezed strongly there. To 
the contrary, in the disk, the gas seems to move across 
the magnetic field line; the gas squeezes and drags the 
magnetic field lines to the center. It is to be noted 
that the magnetic field strengths in and out of the 
disk {p ^ 10^ or p <; 10^) differ slightly, although 
the gas densities differ in three orders of magnitude. 
This is very similar to the flow realized in a collapsing 
rotating isothermal cloud (Norman, Wilson, & Bar- 
ton 1980). That is, the gas accreting from the direc- 
tion of angular momentum moves mainly along the 
lines of constant specific angular momentum; the gas 
contracting radially in the disk seems to move across 
these lines. 

Comparing Figures 4a and &, it is shown that the 
thickness of the disk decreases according to the in- 
crease of Pc'., The half thickness of the disk decreases 
from 5 X 10^'* (a) to 1.5 x 10^^ (&), while the central 
density increases from 10^'^ to 10^°. When pc reaches 
10^", although the resolution obtained by L9 (b) is in- 
sufficient, the Lll grid resolves the vertical structure 
of the disk even in the central high-density region. 

In Figure 5, the cross-cut view along the z-axis (a) 
and that of r-axis (b) are plotted. From Figure 5a, 
it is clear that the cylindrical cloud is divided into a 
disk and an intra-disk. A half-height of the disk de- 
creases as from h ^ 10^^ for pc ^ 10"'' to h ^ 10^'* 
for Pc ^ 10^". Gas accrets to the disk with super- 
sonic speed jw^l > 1. As long as the gas density in 
the disk does not exceed pc ~ 10^, the infalling gas is 
gradually decelerated and forms a disk. However, af- 
ter the disk is developed pc ^ 10^, the gas is stopped 
abruptly on the disk surface where the gas is com- 
pressed and the density increases much. Comparing 
Bz{z, 0) and p(z, 0), in the intra-disk region the mag- 
netic field strength is almost proportional to the den- 
sity, e.g., i?z(z,0) ^ v4p(z,0). The ratio of the ther- 
mal pressure to the magnetic pressure, /3, decreases 
from /3 ~ 1 at z = ±^2/2 to (3 ^ 10^^'^ on the surface 
of the disk (the values are for pc = 10^). This explains 
why the inflow in intra-disk region is controlled by the 
magnetic field and the gas flows along the magnetic 
fields. The j3 value begins to increase after the gas 



enters the disk and /3 reaches again ^ 7.5 reaching 
the center. Thus the magnetic fields are dragged by 
radial motion of the gas. 

To the contrary, there is no discontinuity in the ra- 
dial distributions (5). The density increases in a sys- 
tematic fashion. There seems to exist an asymptotic 
density distribution and the radial density distribu- 
tion reaches it. It is well fitted as yo(0, r) ~ 20r~^-''^. 
The power law distribution ex r~^ is seen in a course 
of collapse of the isothermal rotating cloud (Norman 
et al 1980; Narita et al 1984) as well as the isother- 
mal spherical collapse (Larson 1969; Pension 1969). 
The disk is divided into two parts: a core which 
shows /9(0, r) ^constant and an envelope which shows 
p ex r^^. The evolution proceeds as a fashion that the 
core size Vc decreases monotonically and the region of 
envelope predominates. 

The magnetic field strength seems to approach to 
a singular distribution as i3z(0, r) ~ 2r^^°^. It is 
to be noted that the derived power is nearly equal 
to -1 as Bz{0,r) ex r~^, which is a consequence of 
a tight relationship between /o(0, r) and i?2(0, r) as 
5(0, r) ex p(0, r)i/2 (see Fig.5c; Tomisaka 1995; Naka- 
niura et al 1994). Since the magnetic fields are frozen- 
into the gas, the column density integrated along the 
magnetic field line is proportional to the magnetic 
field strength. If the disk is approximated as an 
isothermal plane- parallel disk, the column density a 
is related to the midplane density po ss, 

ttG 2 /„ 1 \ 

pa = TT^cr +ps, (3-1) 

where ps means the surface density = Pcxt/c^ (Spitzer 
1978). This indicates that po oc er^ when po S> Ps- 
Thus, in the plane-parallel approximation, the mag- 
netic field strength 3^(0, r) is proportional to the 
square root of the density p(0, z). 

The column density which is integrated vertically 
as cr = /_; %2 P'^^ ^^ plotted in Figure 5d. It shows 
that the distribution reaches a singular solution as 
a ex r^^, although it shows a flat distribution in an 
early phase. The singular power law distribution is 



approximated as cr ~ 15r 



From Figure 5d, it is 



shown that the core size, r^ which is deflned as the re- 
gion where a ~ constant, decreases with time and the 
column density in the core CTc is proportional to l/vc, 
that is, (Tc X Tc ~ constant. This is seen in a course 
of "run-away collapse" in rotating isothermal clouds 
(Norman et al 1980; Narita et al 1984). In the run- 
away collapse, the density in a small part of the cloud 



increases infinitely, even the centrifugal force prevents 
the cloud from global contraction. In contrast to 
the angular momentum, the magnetic field can not 
support a cloud as a whole when the mass is larger 
than the critical mass (Mouschovias & Spitzer 1976; 
Tomisaka, Ikeuchi, Nakamura 1988, 1989). Although 
the angular momentum and the magnetic fields have 
different effects on the global contraction, it is to be 
noted that the high density core in the disk cloud con- 
tracts as a 'run-away collapse' in isothermal clouds. 
As a result, it is shown that the radial distributions 
are reaching singular form as 



p{r) ex Bz{r) oc cr{r) ex r 



(3-2) 



The radial inflow velocity is supersonic after p^ > 
10^. However, it is to be noticed that the radial ve- 
locity is much slower than the vertical velocity. This 
is due to the effect of the magnetic fields; that is, 
the magnetic pressure is comparable with the ther- 
mal pressure in the disk and thus the motion across 
the magnetic fields is more or less blocked. 

3.1. Weak Magnetic Fields 

Model B corresponds to a case with weak magnetic 
fields [a = 0.1). Since the magnetic fields support the 
cloud in the r-direction, the radial size of the cloud is 
reduced as i? ~ 0.8694 compared with R ~ 1.039 of 
Model A. The structures for three stages are shown 
in Figure 6. The difference from Model A is appar- 
ent if compared Figure 6a (t = 1.612) with Figure 2a 
{t = 1.367), which has the similar pc- At that epoch, 
the shape of a region which has a positive density 
fluctuation {p > 10^) is a prolate spheroid, which co- 
incides with the eigen-function derived by the linear 
perturbation theory. A disk is not formed yet but a 
spherical core is formed. It has been shown in Paper 
I that in a non-magnetized cylindrical cloud a spher- 
ical core is formed and it continues to collapse. If 
the cloud contracts only in the r-dircction, the mag- 
netic pressure increases in proportional to the square 
of the thermal pressure. Even an initial magnetic 
pressure is as small as 0.05 x thermal pressure, it be- 
comes comparable to the thermal pressure when the 
central density exceeds 10^. Thus before the epoch 
when Pc ^ 10^, the magnetic field only has a sub- 
sidiary role in the dynamics of the cloud. 

In Figure 6b, the next stage is shown. The cen- 
tral density reaches ^ 10^ at that time and the mag- 
netic pressure has an important role compared with 



the previous stage. It is seen that a small disk is 
formed which runs perpendicular to the magnetic 
fields. Comparing this with Figure 36, both of which 
have similar central densities round ^ 10^, the dif- 
ference in the structures of the formed disks is clear. 
It is evident that the lateral size of the disk becomes 
much smaller in case of weak magnetic fields. A sub- 
stantial part of the high-density region still forms a 
sphere even after a small disk contracts and the den- 
sity reaches ~ 10^" (Fig. 6c). 

The cross-cut views along z- and r-axis are shown 
in Figure 7. The respective lines correspond to dif- 
ferent epochs. The characteristic points are similar 
to Model A; magnetic field strength 5^(0, r) is pro- 
portional to the square root of the gas density. It is 
shown that the distributions of density and magnetic 
fields reach asymptotically as p{0, r) ~ 2.5r~^'* and 
Bz{0,r) ~ 0.32r~i-26. The column density distribu- 
tion is well fitted by a power law singular solution as 
a = 24r-0-«74^ 

In Figure 8, the plasma /3, which is the ratio of the 
thermal pressure to magnetic pressure, is plotted ver- 
sus the distance from the center of the disk. Figures 
8 a and b correspond to cross-cuts along the z-axis for 
Model A and B, respectively. With reaching the cen- 
ter (3 decreases and just outside the disk /3 <; 10^^ 
when Pc ^ 10^ for Model A. The minimum of /3 is af- 
fected by the flux contained in the cloud; in Model B 
it is /? ~ 0.1 at an epoch of pc — 10^". In the disk, (3 
increases with reaching the center. In the intra-disk 
region, the thermal and the magnetic pressures are 
well correlated as 91ogp/91ogpinag — 0.6 (for Model 
A) and 0.75 (for Model B). 

Figures 8 c and d correspond to cross-cuts along the 
r-axis. In contrast to /3(z,0), /3(0,r) increases with 
time in the disk. For example, in the range of r ^ 
10~^ the P increases with time. Noting that a = 1 is 
equivalent to /3 — 2, it is clear that the outer part does 
not experience violent contraction. In Model A, the 
j3 seems to reach asymptotically ~ 10. In Model B 
this is found in the range of 15 - 35. Since (3 decreases 
when the magnetic fields are compressed by the radial 
contraction but increases when the cloud contracts 
along the magnetic fields, the steady increase of /3 
indicates that disk formation continues. 

Since the central column density ac{r ~ 0) is pro- 
portional to the magnetic field strength B^ , it is esti- 
mated as 
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where we used equation ( |2-l| ) and the initial condition 
(the suffix "init" represents the initial value). Using 
equation (3-1), the plasma (3 in the central region of 
the cloud is estimated as 



Pc 



507^/3c init, 

1007Va- 



(3-5) 
(3-6) 



This expects fie ^ 100 - 1000 for Models A and B. 
However, the actual Pc ^ 10 is much smaller than 
that expected from a thin plane-parallel disk model 
/3c ~ 100 — 1000. This suggests the gravity gz may be 
overestimated in the plane-parallel disk model. Since 
the disk has a finite radial size, in other word, the col- 
umn density cr(r) decreases proportionally to ex r^^ , 
the gravity should be weaker than that expected by 
a perfectly plane-parallel disk. 

In the isothermal plane-parallel disk, the density 
distribution is expressed as 



P(^) 



sech 



Cs/^/2-kGpc 



(3-7) 



where Pc represents the density on the mid-plane 
z = (Spitzer 1978). Comparing this with p{z,0) 
of Figure 5 a, it is seen that the size is underesti- 
mated if we use the plane-parallel disk approxima- 
tion: When the size is measured by the z-coordinate 
where p{z,0) = 10^, equation (3-7) gives the size of 



~ 2/3 times smaller than the actual size in Figure 5a 
of Pc = 10^°. If we compare the size of p{z, 0) = 10*, 
the factor is equal to ~ 0.4. This indicates that the 
gravity in the plane-parallel disk with an identical 
central density is stronger than that realized in the 
contracting disk here. 

4. Discussion 

4.1. Effect of the Equation of State 

In the previous section, it is shown that the cylin- 
drical cloud forms disks and the central part of the 
disk experiences the "run-away collapse" . However, 
this seems to due to the assumption of the isother- 
mal equation of state for the gas. Actually, before 
the gas density increases infinitely, the cloud core be- 
comes optically thick and the radiation, which keeps 
the gas isothermal, begins to be trapped in the core 
(Hayashi 1966). The critical density, above which the 



equation of state is no more isothermal, is estimated 
Pcrit ~ 10"^"' — 10~^^'^gcm~'^ for the mass range of 
IOOMq — O.OlAfQ. The gas obeys another equation 
state above Pcrit, which should be solved by the radia- 
tive magnetohydrodynamics. Here, for simplicity, we 
assume the polytropic relation between the pressure 
and the density as 



P 



cIp 



ctPcritif/Pcrit) 



a p< Petit, .. ,N 

\I p > Pcrit, 



where A/" means the polytropic index. In Model C 
we assume A/" = 3 and pcrit — lO^ps- If we adopt 



Ps = lOOcm , this means pcrit 



lO^^cni- 



4.1.1. Model C 

Since the difference is expected in the evolution 
after the maximum density is larger than Pcrit, we 
will focus on the central part of the disk. Figure 
9a-c shows the structure of the central disk. It is 
noticed that there exists a steep density jump near 
z ^ 2y. 10""* and the accreting flow in the ^-direction 
is stopped here. The position of density jump seems 
stationary after the central density reaches pc — 10^, 
while the shock front recedes continuously in isother- 
mal models (Compare Fig. 5a with Fig. 9a). In con- 
trast to the cross-cut view along the x-axis, there 
is a slight change in the radial distribution of den- 
sity. However, it is shown that the magnetic fields 
Bz is compressed by the effect of the thermal pres- 
sure, which is relatively higher than Model A (Fig. 5c 
& Fig. 9a). The plasma /3 at the center of the disk 
decreases from ~ 6 of Model A to ~ 3. 

Another difference is the central spherical core with 
a radius of r ~ 1.5 x 10"''. This spherical core keeps 
its size. In the case of isothermal clouds (Models A & 
B), the accumulated mass results in the strong grav- 
itational field as well as the pressure force. To the 
contrary, the pressure force becomes relatively impor- 
tant over the gravity due to the harder equation state 
in this case. Thus, the thermal pressure, which is 
isotropic in nature, plays an important role in the 
core and forms a spherical core. If we assume the 
polytropic index J\f smaller than 3, (that is, a harder 
equation of state) , we may find a bounce due to the 
thermal pressure similar to the cases of the spherical 
cloud (Larson 1969) and the rotating cloud (Narita et 
al. 1986). 



4.1.2. Model D 

The structure in the spherical core with the poly- 
tropic index A/" = 3 is simple for Model C as long as 
Pc <; lOOpcrit- To see the core evolution further, a 
model with pcrit — 10^ is studied. The evolution is 
traced in the range of 4 orders of magnitude of the 
central density, 10^ ^ Pc <, 10^". The polytropic in- 
dex is unchanged. 

Figure 10 shows the structure of the central \z\ <j 
2 X 10-3 (a) and |z| ;$ 4 X IQ-"* (6). Two shock front 
systems are found: one is seen \z\ ^ lO^^ near the 



z-axis. The other is 



1.5 X 10 near the 



axis, which seems to be connected the density jumps 
seen in r ~ 4 x lO""*, \z\ <, 2.5 x lO""*. The first 
shock is related to the density jump appeared in the 
previous model. The structures for pc — 10^, 10^, 
..., 10^" are shown in Figure 10c and d. It shows 
that a jump in Vz appears after pc > 10^. At first, 
the infalling velocity v^ is completely stopped through 
the shock front and the core seems static {pc ^ 10^). 
However, after pc ^ 10* , the core begins to contract 
and the inflowing velocity in the core is accelerated. 
Finally, another density/velocity jump appears when 
pc ~ 10^°. The second structure is seen in Figure 106 
as a region where the density gradient is steep. 

Oblique shocks are seen near r ~ 4 x 10^^. Mat- 
ter flowing from upper-right and upper-left directions 
{\vz\ ~ \vr\) changes its direction through this sur- 
face to the radial direction {\vr\ ^ \zz\). Plotting the 
magnetic fields together, it changes its direction as 
the same sense as the velocity field. 

It is concluded that the collapse in the center of the 
cloud/disk continues even after the equation of state 
becomes effectively harder than the isothermal one. 
As the contraction proceeds, the infall velocity is ac- 
celerated and the second discontinuity is formed. Dis- 
continuity facing to the radial direction is also formed. 
As a course of the non-isothermal contraction, spheri- 
cal or quasi-spherical core is formed. This shows that 
the smooth "run-away" collapse seems to be seen in 
the isothermal collapse and in the actual interstel- 
lar environment, shocks are inevitably formed in the 
dense core. 

4.2. Evolution of a Cylindrical Cloud 

It is shown that the cylindrical cloud threaded by 
magnetic fields running in the direction parallel to the 
symmetric axis is fragmented by the gravitational in- 
stability (Paper I and the present paper). The frag- 



ment which takes initially prolate spheroidal shape 
contracts along the magnetic field and forms a disk. 
Finally, the central part of the disk begins "run-away 
collapse" . Is this a common evolution track? 

To answer the question, we have to start with the 
result of the linear instability theory. Nakamura et 
al. (1993) obtained eigen-functions and growth rates 
of the perturbation added to the isothermal cylinder. 
The wave-length which has the maximum growth rate 
is approximated as ~ 20^ Cs / ^/4^TGp^. The fact that 
the expression does not contain the magnetic field 
strength is important. This means that even when 
the cylindrical cloud is supported in the radial direc- 
tion by the magnetic field, the cloud fragments with a 
characteristic size of scale-length determined by only 
the isothermal sound speed, Cg, and the central den- 
sity, Pc. The fragmentation process is driven by the 
Jean instability and the Jean mode is activated by 
the motion in the z-direction. However, the magnetic 
fields have no effect to the motion parallel to it. This 
is the reason why the maximum growth wave-length 
is not related to the magnetic fields. (Precisely speak- 
ing, the maximum growth wave-length is dependent 
on a as Amax oc a^'^ for a 3> 1. However its depen- 
dence is weak [see Nakamura ct al 1993].) 

As shown in Paper I, if the fragmented disk has 
enough mass-to-flux ratio, M/<I> = cr/i?z, the re- 
sultant disk is supercritical and must contract as a 
whole. The critical value of the mass-to-fiux ratio 
is M/($/G'i/2) ~ 0.17 ~ l/(27r) (Mouschovias & 
Spitzer 1976; Tomisaka, et al 1988, 1989). Since the 
mass-to-flux ratio for the disk formed in the cylindri- 
cal cloud is estimated as 
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(4-2) 



(equation[4.4] of Paper I), models that we calculated 
in the present paper lead to supercritical disks. It is 
reasonable that the "run-away" collapse is seen in all 
models. 

How about a cloud with strong magnetic fields? 
Since a increases as the cloud contracts radially, it 
is realized as a consequence of contraction. When 
hydrostatic balance is achieved after contraction, the 
gravitational instability grows. If the value of a at 
that time is over ^ 100, the mass-to-fiux ratio of 
the first growing fragment may be smaller than the 
critical value. However, it should be noted that the 
fragment which is more massive than the most unsta- 
ble perturbation can grow even with a smaller growth 



rate than the most unstable one. This means the frag- 
ments tend to merge each other moving parallelly to 
the magnetic fields. This makes the mass increase but 
keeps the magnetic flux unchanged. As a result, even 
when the cloud is compressed and a in the central 
part of the cloud is over 100, it is unlikely that many 
disks are formed and exist stably. 

4.3. Application to Dimensional Values 

The physical values were scaled with units of c^ ~ 
47rG = Pcxt = 1- In this subsection, let us convert 
the non-dimensional values to dimensional ones. As- 
suming the temperature of T = lOK and molecular 
weight of /i = 2.33, the isothermal sound speed for 
this temperature is Cs = 190ms-i(r/10K)i/2. This 
leads to length scale of 
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,100cm-3y V200ms-i/^ " -* 

where Ug means the particle density on the surface 
of the cloud. Since the radii of the initial states of 
Models A and B are equal to R = 0.8694 and 1.039 
respectively, the radial size of the cloud corresponds 
to - 0.4pc(ns/100cm-3)-i/2(c,/200ms-i). The ini- 
tial cloud has the density contrast of _F = Pc/ps- 
Typical densities in the cloud center varies from 
Uc = 10'*cm"3(-^/j^oO)(n^/100cm"3) ^^ lO^^cm^s 

(pc/lO^'^Ps) (ris/lOOcm^"^). The magnetic field is nor- 
malized by a unit of 
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V4^rp^, 
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1/2 



VlOOcm-3 
Time is measured with a unit of 



1/^4ttGps, 
1.75 Myr 



200ms- 



-1/2 



- (4-4) 
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Since perturbations in the cylindrical cloud grow and 
form disks in a time scale of t ~ 1.3 — 1.6 for Models 
A and B, the physical time scale for the cylindrical 
cloud to experience the run-away collapse is estimated 
as ~ (2 - 3)Myr(ns/100cm-3)-i/2. 

The separation between the disks is estimated 



0.72 pc 7 
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This indicates that the disk separation is approxi- 
mately equal to the diameter of the cylindrical cloud. 
Since the line density of the cloud is equal to 



A = 



2'Kprdr 



(2c2/G)(F1/2_i)/fV2(i + c./2), 
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for F = 100 and a = 1, the mass of the resultant disk 
is estimated as 



M = l,X, 



18 M07 



/ ns 



200m s-i/ \100cm-3 



-1/2 

(4.8) 



Singular power law distributions of the density {p{Q, r) • 
20r-208), the magnetic field {B^ ~ 2r-i°^), and the 
column density {a ~ 15r~^°^) for Model A are writ- 
ten as follows: 



p{r) ~ 2 X 10^ cm 
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Recently, it is found that a gas disk around HL 
Tauri is contracting towards the center (Hayashi, 
Ohashi, & Miyama 1993). The disk with a mass of 
Afdisk ~ 0.022 - O.IIMq in i? = 1400 AU shows an 
indication of inflow with a speed of u^ ^ 1 kms~ . 
If this object is a disk which is in a stage of the 
run-away coUase we studied, the observed contrac- 
tion speed inevitablly leads to a large sound speed 
Cs ~ 0.5km s~ , because the radial inflow speed is 
approximately equal to ~ 2cs- But, since this gives 
a relatively large disk mass as iWdisk ~ /g (J^irrdr ~ 
lMQ(c^/0.5kms"^)(nJ100cm-3)i/2(ij/o.36pc)i02in 
R =1400 AU, HL Tauri seems an object in a free-fall 
stage or the inside-out collapse stage (Galli & Shu 
1993a, b) subsequent to the run-away collapse. Thus, 



if objects are found with relatively low contraction 
speed Vr ~ 2cs, this might be an object just before 
the run-away collapse. 

This work was supported in part by the Grants- 
in-Aid from the Ministry of Education Culture and 
Science (05217208) in fiscal years 1993-1994. Numer- 
ical calculations were done in part by supercomputer: 
S-3800/480 in University of Tokyo and VPP500/7 in 
IS AS. 

A. Nested Grid Scheme 

In Appendix, the nested grid scheme is described 
very briefly. We used 15 levels of the grids, in which 
finer grids are prepared for the very center of the con- 
tracting disk and coarser grids are to calculate the 
whole cloud. As seen in Figure 1, a finer grid is wholly 
covered by a coarser grid. We name respective grids 
as LO, LI, . . ., L14. LO is the coarsest grid and L14 
is the finest. In the scheme, we choose a mesh size 
of a finer grid (Ln) as 1/2 of that of the coarser grid 
(Ln — 1). Since the mesh number of the respective 
grids are taken the same (64) , the size of the numer- 
ical box of the coarser grid (Ln — 1) is twice as large 
as that of the finer one (Ln) . For the character of the 
problem, all multiple grids are put in a co-centered 
fashion; the coordinates of the centers of the z-axis 
are taken identical for all grids (Fig.l). 

Basic equations are the unste ady MHD equation 



and the Poisson equation (eq. 2-5 - 2-10 ). As a 



boundary condition for Ln, quantities on the upper 
{z = Zu), lower (z ~ Zi), and outer (r — Ro) bound- 
aries, such as p, V, B, and ip, are calculated from 
those of Ln — 1. To derive the boundary values, we 
interpolate from those on coarser grids with van Leer's 
monotonic interpolation (van Leer 1977). After deter- 
mining the boundary values, we can solve the MHD 
and Poisson equations for Ln. Then, equations for 
Ln — 1 are solved. Since one quarter of the mesh 
points of Ln — 1 are covered by the Ln grids, quanti- 
ties of Ln — 1 in the overlapped region are overwritten 
by new values of Ln. To derive quantities of Ln — 1 
we simply take the volume-average of the quantities in 
corresponding 4 cells of Ln. Quantities in the region 
overlapped by Ln — 1 are calculated but are not actu- 
ally used, when equations for Ln — 1 are solved. This 
makes data structure simple and the code achieves 
high performances on vector machines. 

We show the algorithm here for three grids: coarse 



(LO), middle (LI) and a fine grid (L2), for simplicity. 
(0.0) Calculate boundary values for LI from LO 

(1.0) Calculate boundary values for L2 from LI 

(2.1) Time step calculated from Courant 

condition {dt21) 
f2.2^ Solve MHD eq. for L2 
(2.3) Solve the Poisson eq. for L2 

(2.1) Time step calculated from Courant 

condition {dt22) 
f2.2) Solve MHD eq. for L2 
(2.3) Solve the Poisson eq. for L2 

(1.1) Time step for LI is set as {dtll = dt21 + 
dt22) 

(1.2) Solve MHD eq. for LI 

(1.3) Solve the Poisson eq. for LI 

(1.4) A part of LI overlaid by L2 are replaced 
by quantities of L2 

(1.0) Calculate boundary values for L2 from LI 

(2.1) Time step calculated from Courant 

condition {dt2\) 
f2.2^ Solve MHD eq. for L2 
(2.3) Solve the Poisson eq. for L2 

(2.1) Time step calculated from Courant 
condition {dt22) 

2.2) Solve MHD eq. for L2 

2.3) Solve the Poisson eq. for L2 

(1.1) Time step for LI is set as {dtl2 = dt21 + 
dt22) 

(1.2) Solve MHD eq. for LI 

(1.3) Solve the Poisson eq. for LI 

(1.4) A part of LO overlaid by LI are replaced 
by quantities of LI 

(0.1) Time step is set as dtQ = dtll + dtl2 
(0.2) Solve MHD eq. for LO 
(0.3) Solve the Poisson eq. for LO 
(0.4) A part of LO overlaid by LI are replaced by 
quantities of LI 

As is seen, in the n-th level, we calculate the MHD 
and Poisson equation of Ln + 1 twice and then go to 
Ln. The time step of Ln is equal to the total of two 
time steps of Ln -1-1. In the present paper, the time 
step is determined by the finest grid LiV and used in 
coarser grids LO, LI, . . ., LiV— 1. The finest grid level 
are changed with necessity; that is, at first N is equal 
to 4; and a finer grid level is generated {N — > TV -I- 1) 
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when a density contrast in the finest grid becomes 
larger than 10. 

We solved the Poisson equation by the "modified 
incomplete LU decomposition preconditioned conju- 
gate gradient squared method (MILUCGS: Meijerink 
& van der Vorst 1977; Gustafsson 1978). For MHD 
equations we use "monotonic scheme" (van Leer 1977, 
Norman & Winkler 1986) and "constrained transport 
scheme" (Evans & Hawley 1988). Routines to solve 
MHD equations (0.2, 1.2, & 2.2) in the above are ac- 
tually an identical subroutine as well as the Poisson 
solver (0.3, 1.3, & 2.3) are identical with one another. 
Since the mesh numbers of each levels are taken as 
the same, we can use the identical routine for data 
belonging to different levels of grids. Due to this sim- 
plicity, the finest grid level (N) is changed easily with- 
out changing the code. The code for a single level of 
grids is the same as the previous paper (Tomisaka 
1995). 
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Fig.l. — Explanation of the nested grid method. Level covers the entire cloud. Level 1 covers the central 
1/4 of LO, Level 2 covers the central 1/4 of LI, and so on. Since the numbers of the meshes of respective levels are 
the same (64x64), the LI has twice higher resolution than LO. 

Fig. 2. — Snapshots of Model A at the epoch of pc = 10'"' [t = 1.367). The Left panel (a) shows the structure in 
LO (the coarsest grid) and the right panel (&) shows that of L2 (four times closer than LO). Density contour lines 
(thick lines) and magnetic field lines (thin lines) are illustrated. The contour levels of the denisties are shown by 
numbers attached to the contour lines as logj^Q p. Numbers seen on the right-bottom corner represent the maximum 
and minimum values of logp appeared in respective panels. 

Fig. 3. — Two snapshots of Model A: the left panel corresponds to the epoch of Pc — 10® {t — 1.372) (a) and 
the right shows that of pc = 10^'^ {t = 1.373) (b). The fifth level of grids L5 are plotted. Density contours and 
velocity fields are illustrated. Velocity vectors are illustrated every 3x3 meshes. Numbers seen above the density 
range means the maximum value of |v| (3.36 left and 4.02 right). 

Fig. 4. — The upper panels: L9 structures at the epochs of Pc = lO*'^'' (left: a) and pc = 10^"-^^ (right: b) for 
Model A. The lower shows the state when p^ — 10^" but in the Lll grid (c). Density contour lines (thick lines), 
magnetic field lines (thin lines), and velocity fields (vectors) are illustrated. 

Fig. 4. — (a): Cross-cut view on the z-axis for Model A. p{z,0) (solid lines), Bz{z,0) (dashed lines), and 
|uz(2;,0)| (dotted lines) are illustrated. Eight different epochs are snapshoted. pc — lO'^"^, 10**, lO'''"^, 10®°^, 
lO"^-^, 10^-^^, 10^-^2 and lO^o^^ are plotted. (6): Cross-cut view on the r-axis. p{0,r) (solid lines), B^(0, r) (dashed 
lines), and —Vr{0,r) (dotted lines) are shown. The epochs of the snapshots are the same as a. (c): The correlation 
between density and magnetic field strength on the disk midplane z = 0. Eight lines of /9(0,r)-i?2(0,r) correspond 
to the respective epochs of a. The left-bottum corner is the surface of the cloud and the right-upper ends are the 
disk centers, (d): The column density distribution against r. cr(r) = J pdz. The epochs of the snapshots are the 
same as a. 

Fig. 6. — The structure of the fragment in the cylindrical cloud for Model B. (a): the snapshot at i = 1.612 
{pc — 10^'^^) in L2. High-density part of the cloud indicates a prolate spheroidal shape and it contains a spherical 
core. (6): the snapshot at t = 1.616 {pc = 10^'^^) in L5. A small disk is formed in the center of the spherical core, 
(c): the snapshot at t = 1.616 {pc = 10i°°i) in Lll. 

Fig. 7. — The same as Fig. 5 but for Model B. Eight different epochs are snapshoted: pc ~ 10'^°'^, 10^'°^, 10^-^^, 
106-5, 10717^ 108-38, 109-19, and 1010-01. 

Fig. 8. — The plasma /3 versus the distance from the center of the disk \z — Zc\ for Models A (a) & B (6). 
Respective curves correspond to different epochs. The epochs of the snapshots are the same as Figs. 5 & 7. As the 
contraction proceeds the minimum of the /3 decreases. The plasma /3 versus radial distance r for Models A (c) & B 
(d). Generally, the (3 in the disk increases with time. It indicates that the mass accretion continues on to the disk. 

Fig. 9. — Model C, in which the equation of state is assumed hard for p > lO^p^. This mimics the situation that 
the gravitational energy liberated in a course of collapse cannot be transferred after the cloud becomes optically 
thick, (a): snapshot when pc = lOioo^. Density contour lines, magnetic fields, and velocity fields are plotted. (6) 
and (c): the cross-cut views along the z- and r-axis, when pc = 10^-^ (i = 1.373), lO^i^ (t = 1.374), and lOiO-oe 
(t= 1.374). 

Fig. 10. — (a) & (6): snapshots of Model D (t = 1.374 & Pc = lO^o^^), which is similar to model C but the 
critical density above which the equation of state is hard is assumed lO^Ps- Comparing level 9 (a) and level 11 (6), 
it is shown that two shock fronts exist: the outer one is seen in (a) in a distance of r ~ 10"^ from the center of 
the cloud core. Another one is seen in (6) at z ~ 1.5 x 10""* near the z-axis and r ~ 3.5 x 10^^ and z ~^ 0. (c) 
& (d): cross-cut views along the z- and r-axis, respectively. Snapshots correspond to pc = 10® (t = 1.372), lO^i 
(t = 1.373), 108-21 {t = 1.374), 109, {t = 1.374) & lOioos (t = 1.374) 
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